 $ontext
[1]        $model:sj
[2]        
[3]        * Stylized Johansen model
[4]        
[5]        $sectors:
[6]          xcom(SECT)              ! production
[7]          y                 ! aggregate consumption
[8]        
[9]        $commodities:
[10]         pc(SECT)              ! price of commodity
[11]         py                ! price of aggregate consumption
[12]         pf(FAC)             ! price of primary factor
[13]       
[14]       $consumers:
[15]         ra                ! representative agent income
[16]       
[17]       $prod:xcom(SECT)
[18]       * Tricky. Can infer that SECT1 is a subset of SECT here since pc is defined over SECT
[19]       *  but is also used over SECT1.
[20]         o:pc(SECT)         q:DVCOM(SECT)
[21]         i:pc(SECT1)        q:DVCOMIN(SECT1,SECT)
[22]         i:pf(FAC)          q:DVFACIN(FAC,SECT)
[23]       
[24]       $prod:y
[25]         o:py                 q:200
[26]       * NB. This error not picked up ** Need to do so in Tom's part of code??
[27]       * Error:  i:p(SECT)            q:DVHOUS(SECT)
[28]         i:pc(SECT)            q:DVHOUS(SECT)
[29]       
[30]       $demand:ra
[31]         d:py
[32]         e:pf(FAC)           q:endow(FAC)
[33]       
 $offtext

*   Alert the user to undefined sets or parameters which are referenced in the MPSGE model:

$if not settype sect        $abort "sect should be a SET."
$if not settype fac         $abort "fac should be a SET."
$if not settype sect1       $abort "sect1 should be a SET."
$if not partype DVCOM       $abort "DVCOM should be a PARAMETER."
$if not dimension 1 DVCOM   $abort "DVCOM has an inconsistent dimension."
$if not partype DVCOMIN     $abort "DVCOMIN should be a PARAMETER."
$if not dimension 2 DVCOMIN $abort "DVCOMIN has an inconsistent dimension."
$if not partype DVFACIN     $abort "DVFACIN should be a PARAMETER."
$if not dimension 2 DVFACIN $abort "DVFACIN has an inconsistent dimension."
$if not partype DVHOUS      $abort "DVHOUS should be a PARAMETER."
$if not dimension 1 DVHOUS  $abort "DVHOUS has an inconsistent dimension."
$if not partype endow       $abort "endow should be a PARAMETER."
$if not dimension 1 endow   $abort "endow has an inconsistent dimension."

*    Define aliases for sets:

alias (sect,sect_);   
alias (fac,fac_);     
alias (sect1,sect1_); 

*    Benchmark costs and value shares:

PARAMETER

*    Benchmark costs:

        c0_002(sect,sect1) Benchmark cost for input pc(SECT1) in sector xcom(SECT)
        c0_003(sect,fac)   Benchmark cost for input pf(FAC) in sector xcom(SECT)
        c0_001(sect)       Benchmark cost for input nest s in sector xcom(SECT)
        c0_007(sect)       Benchmark cost for input pc(SECT) in sector y
        c0_006             Benchmark cost for input nest s in sector y

*    Benchmark revenue:

        r0_005(sect) Benchmark revenue for output pc(SECT) in sector xcom(SECT)
        r0_004(sect) Benchmark revenue for output nest t in sector xcom(SECT)
        r0_009       Benchmark revenue for output py in sector y
        r0_008       Benchmark revenue for output nest t in sector y

*    Benchmark shares:

        theta_002(sect,sect1) Benchmark share for input  pc(SECT1) in nest s for sector xcom(SECT)
        theta_003(sect,fac)   Benchmark share for input  pf(FAC) in nest s for sector xcom(SECT)
        theta_005(sect)       Benchmark share for output pc(SECT) in nest t for sector xcom(SECT)
        theta_007(sect)       Benchmark share for input  pc(SECT) in nest s for sector y
        theta_009             Benchmark share for output py in nest t for sector y
;



*   Invoke the function preprocessor


$sysinclude gams-f

c0_002(sect,sect1) = DVCOMIN(SECT1,SECT);                                
c0_003(sect,fac) = DVFACIN(FAC,SECT);                                    
r0_005(sect) = DVCOM(SECT);                                              
c0_001(sect) = +sum(sect1,c0_002(sect,sect1))+sum(fac,c0_003(sect,fac)); 
r0_004(sect) = +r0_005(sect);                                            
c0_007(sect) = DVHOUS(SECT);                                             
r0_009 = 200;                                                            
c0_006 = +sum(sect,c0_007(sect));                                        
r0_008 = +r0_009;                                                        
theta_002(sect,sect1) = c0_002(sect,sect1)/c0_001(sect); 
theta_003(sect,fac) = c0_003(sect,fac)/c0_001(sect);     
theta_005(sect) = r0_005(sect)/r0_004(sect);             
theta_007(sect) = c0_007(sect)/c0_006;                   
theta_009 = r0_009/r0_008;                               

*   Cost functions

p_002(sect,sect1) == pc(SECT1);                                                                                     
p_003(sect,fac) == pf(FAC);                                                                                         
p_005(sect) == pc(SECT);                                                                                            
p_001(sect) == sum(sect1,theta_002(sect,sect1) * p_002(sect,sect1))+sum(fac,theta_003(sect,fac) * p_003(sect,fac)); 
p_004(sect) == theta_005(sect)*p_005(sect);                                                                         
p_007(sect) == pc(SECT);                                                                                            
p_009 == py;                                                                                                        
p_006 == sum(sect,theta_007(sect) * p_007(sect));                                                                   
p_008 == theta_009*p_009;                                                                                           

*   Quantity functions

q_002(sect,sect1) == DVCOMIN(SECT1,SECT); 
q_003(sect,fac) == DVFACIN(FAC,SECT);     
q_005(sect) == DVCOM(SECT);               
q_007(sect) == DVHOUS(SECT);              
q_009 == 200;                             

*   Compensated demand indices:

d_001(sect) ==  1;                                                    
d_002(sect,sect1) == d_001(sect)* (p_001(sect)/p_002(sect,sect1))**0; 
d_003(sect,fac) == d_001(sect)* (p_001(sect)/p_003(sect,fac))**0;     
d_006 ==  1;                                                          
d_007(sect) == d_006* (p_006/p_007(sect))**0;                         

*   Compensated supply indices:

s_004(sect) ==  1;                                        
s_005(sect) == s_004(sect)* (p_004(sect)/p_005(sect))**0; 
s_008 ==  1;                                              
s_009 == s_008* (p_008/p_009)**0;                         

*   Endowments:

e_010(fac) == endow(FAC); 


POSITIVE VARIABLES
        xcom(SECT) production
        y          aggregate consumption
        pc(SECT)   price of commodity
        py         price of aggregate consumption
        pf(FAC)    price of primary factor
        ra         representative agent income;


EQUATIONS
        prf_xcom(SECT) "production"
        prf_y          "aggregate consumption"
        mkt_pc(SECT)   "price of commodity"
        mkt_py         "price of aggregate consumption"
        mkt_pf(FAC)    "price of primary factor"
        inc_ra         "representative agent income";



MODEL sj / prf_xcom.xcom,prf_y.y,mkt_pc.pc,mkt_py.py,mkt_pf.pf,inc_ra.ra/; 

*     Level value assignments
xcom.l(SECT) = 1; 
y.l = 1;          
pc.l(SECT) = 1;   
py.l = 1;         
pf.l(FAC) = 1;    

*     Zero profit conditions
prf_xcom(SECT).. c0_001(sect)*p_001(sect) =G= r0_004(sect)*p_004(sect);
prf_y..          c0_006*p_006 =G= r0_008*p_008;

*     Market clearance conditions
mkt_pc(sect).. 
               s_005(sect)*xcom(sect)*q_005(sect)+sum((sect1,sect_)$(sameas(sect,sect1)),-d_002(sect,sect1)*xcom(sect_)*q_002(sect,sect1))-d_007(sect)*y*q_007(sect) =e= 0;
mkt_py..       
               s_009*y*q_009-ra/py =e= 0;
mkt_pf(fac)..  
               +sum(sect,-d_003(sect,fac)*xcom(sect)*q_003(sect,fac)) + e_010(fac) =e= 0;

*     Income definitions
inc_ra.. 
         ra =e=+sum((fac,pf(FAC)*endow(FAC));

*     Auxiliary constraints:
